**A) Test if random slopes are needed for parenting quality
**B) Capture random slopes

global demo "male other_gender secondary tertiary hhinc  age_group2 age_group3 age_group4 age_group5 age_group6 age_group7 foreign_born  employed never_married married live_alone no_children one_child two_children  rural village suburb"
global parent_demo "parents_married living_structure2 living_structure3 living_structure4  child_poverty was_abused"

use "${clean_data}\gfs_cleaned_merged_database.dta", clear
cd "${output}\"

estimates clear

foreach Y in HOPE HAPPY  {
eststo: mixed `Y'  PCRQ $demo $parent_demo  [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
estimates store mod1`Y'
predict rslope_m1`Y' , reffects  relevel(iso)

eststo:  mixed `Y'  PCRQ $demo $parent_demo  [pw=weight] || iso: PCRQ, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
estimates store mod2`Y'
lrtest  mod1`Y' mod2`Y', force
predict rslope_m2`Y' rint_m2`Y', reffects relevel(iso)
}

esttab using "SUP_random_slope_model.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace
 
desc rslope_m1HOPE
desc rslope_m2HOPE
desc rint_m2HOPE

encode country_name, gen(country_num)

gen B1=.
gen B2=.
gen SE1=.
gen SE2=.

forval i=1/22 {
reg HOPE PCRQ $demo living_structure2 living_structure3 living_structure4 [aw=weight] if country_num==`i'
replace B1=_b[PCRQ] if country_num==`i'
replace SE1=_se[PCRQ] if country_num==`i'
}


collapse (mean)  rslope_m2HOPE  rslope_m2HAPPY rslope_m1HOPE  rslope_m1HAPPY rint_m2HOPE rint_m2HAPPY  B1 SE1 v2x_polyarchy zwvs_traditional zwps_index lngdp, by(country)
save "random_slopes.dta", replace

use "random_slopes.dta", clear
cor B1 rslope_m2HOPE  rslope_m2HAPPY rslope_m1HOPE  rslope_m1HAPPY rint_m2HOPE rint_m2HAPPY

reg B1 rslope_m2HOPE rint_m2HOPE
gen c=_b[_cons]
di c
